Dopasuj model (lub modele) przewidujący wartość zmiennej \(ConcreteStrength\) (Wytrzymałość na ściskanie betonu) . Model zbuduj w wybrany przez siebie sposób używając regresji wielorakiej. Oceń jakość modelu i spełnienie założeń. Porównaj otrzymane modele.
Chociaż nie jesteśmy ekspertami w dziedzinie cementu i betonu, zamierzamy zastosować model regresji wielorakiej, by analizować i interpretować dane zawarte w załączonym pliku. Naszym celem jest zrozumienie, jak różne składniki wpływają na wytrzymałość betonu, wykorzystując dostępne dane.
Zbiór danych “Concrete Compressive Strength” to zestaw informacji dotyczących składników i wytrzymałości betonu. Zawiera zmienne odnoszące się do składników mieszanki betonowej oraz ich wpływu na wytrzymałość betonu. Zawiera 9 zmiennych:
Cement (component 1)(kg in a m3 mixture): Ilość cementu w kilogramach na metr sześcienny mieszanki.
Blast Furnace Slag (component 2)(kg in a m3 mixture): Ilość żużla wielkopiecowego w kilogramach na metr sześcienny mieszanki.
Fly Ash (component 3)(kg in a m3 mixture): Ilość popiołu lotnego w kilogramach na metr sześcienny mieszanki.
Water (component 4)(kg in a m3 mixture): Ilość wody w kilogramach na metr sześcienny mieszanki.
Superplasticizer (component 5)(kg in a m3 mixture): Ilość superplastyfikatora w kilogramach na metr sześcienny mieszanki.
Coarse Aggregate (component 6)(kg in a m3 mixture): Ilość kruszywa grubego w kilogramach na metr sześcienny mieszanki.
Fine Aggregate (component 7)(kg in a m3 mixture): Ilość kruszywa drobnego w kilogramach na metr sześcienny mieszanki.
Age (day): Wiek mieszanki betonowej w dniach.
Concrete compressive strength(MPa, megapascals): Wytrzymałość na ściskanie betonu w megapaskalach.
concrete_data <- read.csv("ConcreteData.csv", header=TRUE)
head(concrete_data)
Zmiana nazw zmiennych (na krótsze):
names(concrete_data) <- c("Cement", "BlastFurnaceSlag", "FlyAsh", "Water", "Superplasticizer", "CoarseAggregate", "FineAggregate", "Age", "ConcreteStrength")
head(concrete_data)
Następnie postaramy się dopasować model (lub modele), które bedą przewidywały w sposób liniowy wartość zmiennej \(ConcreteStrength\). Aby to zrobić, potrzebujemy znaleźć zmienną objaśniające, które najbardziej odpowiadają do modelu regresji liniowej wielorakiej. Najpierw zapoznajmy się z naszymi danymi:
Ogólne statystyki:
summary(concrete_data)
## Cement BlastFurnaceSlag FlyAsh Water
## Min. :102.0 Min. : 0.00 Min. : 0.00 Min. :121.8
## 1st Qu.:190.7 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.:166.6
## Median :255.0 Median : 19.00 Median : 0.00 Median :185.0
## Mean :273.9 Mean : 71.07 Mean : 57.87 Mean :181.8
## 3rd Qu.:339.0 3rd Qu.:140.32 3rd Qu.:118.60 3rd Qu.:192.0
## Max. :540.0 Max. :316.10 Max. :200.10 Max. :247.0
## Superplasticizer CoarseAggregate FineAggregate Age
## Min. : 0.000 Min. : 801.0 Min. :594.0 Min. : 1.00
## 1st Qu.: 0.000 1st Qu.: 932.0 1st Qu.:744.3 1st Qu.: 7.00
## Median : 6.700 Median : 968.0 Median :780.0 Median : 28.00
## Mean : 6.141 Mean : 974.4 Mean :777.1 Mean : 38.47
## 3rd Qu.:10.100 3rd Qu.:1030.2 3rd Qu.:822.6 3rd Qu.: 28.00
## Max. :32.200 Max. :1145.0 Max. :992.6 Max. :365.00
## ConcreteStrength
## Min. : 6.27
## 1st Qu.:23.20
## Median :33.42
## Mean :34.98
## 3rd Qu.:44.66
## Max. :82.60
Dokładniejsze statystyki:
summary_data<-skim(concrete_data)
summary_data
| Name | concrete_data |
| Number of rows | 908 |
| Number of columns | 9 |
| _______________________ | |
| Column type frequency: | |
| numeric | 9 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Cement | 0 | 1 | 273.93 | 99.23 | 102.00 | 190.70 | 255.00 | 339.00 | 540.0 | ▆▇▆▃▁ |
| BlastFurnaceSlag | 0 | 1 | 71.07 | 84.65 | 0.00 | 0.00 | 19.00 | 140.32 | 316.1 | ▇▂▂▁▁ |
| FlyAsh | 0 | 1 | 57.87 | 64.58 | 0.00 | 0.00 | 0.00 | 118.60 | 200.1 | ▇▁▃▂▁ |
| Water | 0 | 1 | 181.84 | 19.73 | 121.80 | 166.60 | 185.00 | 192.00 | 247.0 | ▁▅▇▂▁ |
| Superplasticizer | 0 | 1 | 6.14 | 5.50 | 0.00 | 0.00 | 6.70 | 10.10 | 32.2 | ▇▇▁▁▁ |
| CoarseAggregate | 0 | 1 | 974.41 | 76.56 | 801.00 | 932.00 | 968.00 | 1030.25 | 1145.0 | ▂▅▇▅▂ |
| FineAggregate | 0 | 1 | 777.12 | 74.02 | 594.00 | 744.27 | 780.05 | 822.65 | 992.6 | ▂▃▇▃▁ |
| Age | 0 | 1 | 38.47 | 47.54 | 1.00 | 7.00 | 28.00 | 28.00 | 365.0 | ▇▁▁▁▁ |
| ConcreteStrength | 0 | 1 | 34.98 | 16.14 | 6.27 | 23.20 | 33.42 | 44.66 | 82.6 | ▅▇▆▃▁ |
Histogramy dla poszczególnej zmiennej:
par(mfrow = c(3,4))
for (i in 1:ncol(concrete_data)){
column_name<-colnames(concrete_data)[i]
hist(concrete_data[[column_name]],
main = paste("Histogram - ", column_name),
xlab = "Value", ylab = "Frequecy")
}
Możemy przystąpić do wyznaczania zmiennych objaśniających do naszych modeli. Obliczymy współczynniki korelacji Pearsona \(ConcreteStrength\) z każdą poszczególną zmienną w naszej bazie danych oraz stworzymy wykresy pokazujące każdą zmienną do zmiennej \(ConcreteStrength\) (analiza wykresów i korelacji):
correlation_concrete <- cor(concrete_data$ConcreteStrength, concrete_data)
print(correlation_concrete)
## Cement BlastFurnaceSlag FlyAsh Water Superplasticizer
## [1,] 0.4473283 0.1849669 -0.09264568 -0.2786437 0.3595313
## CoarseAggregate FineAggregate Age ConcreteStrength
## [1,] -0.2056304 -0.1488172 0.3969586 1
for (col_name in names(concrete_data)) {
p <- ggplot(data = concrete_data, aes_string(x = col_name, y = "ConcreteStrength")) +
geom_point() + labs(x = col_name, y = "ConcreteStrength") + ggtitle(paste("ConcreteStrength względem", col_name," (Korelacja:", cor(concrete_data$ConcreteStrength, concrete_data[[col_name]]), ")"))
print(p)
}
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Również stworzymy macierz korelacji, aby lepiej również zrozumieć jak wszystkie dane są ze sobą skorelowane:
ggcorrplot(cor(concrete_data), type='upper')
corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)
Najpierw stwórzmy pierwszy model, gdzie będą uwzględnione wszystkie zmienne:
model1 <- lm(ConcreteStrength ~ ., data = concrete_data)
summary(model1)
##
## Call:
## lm(formula = ConcreteStrength ~ ., data = concrete_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.050 -5.071 0.829 5.674 26.235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.627526 26.526229 0.401 0.6888
## Cement 0.116777 0.008190 14.259 < 2e-16 ***
## BlastFurnaceSlag 0.097201 0.009827 9.892 < 2e-16 ***
## FlyAsh 0.075877 0.012028 6.309 4.41e-10 ***
## Water -0.212400 0.041257 -5.148 3.23e-07 ***
## Superplasticizer 0.244278 0.096365 2.535 0.0114 *
## CoarseAggregate 0.005622 0.009247 0.608 0.5434
## FineAggregate 0.007508 0.010535 0.713 0.4762
## Age 0.178619 0.006530 27.352 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.059 on 899 degrees of freedom
## Multiple R-squared: 0.6878, Adjusted R-squared: 0.685
## F-statistic: 247.6 on 8 and 899 DF, p-value: < 2.2e-16
\(R^2\) na poziomie \(0.6878\), oznacza iż model wyjaśnia około 68,8% obserwacji wartości zmiennej \(ConcreteStrength\).
Istnieją wartości p-value blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Największy wpływ na model mają \(Cement\), \(BlastFurnaceSlag\), \(FlyAsh\), \(Water\), \(Superplasticizer\) oraz \(Age\).
Wartość \(p-value\) dla \(Intercept\) większa niż \(0.05\) oznacza, że przecięcie nie jest istotnie różne od zera. Co by onzaczało, że siła wytrzymałości betonu nie potrzebuje żadnych składników. Możliwe zmiany później.
Założenie 1. Liniowa zależność między zmienną objaśnianą, a objaśniającymi:
Możemy stwierdzić na podstawie \(Adj. R^2=0.685\) oraz \(p-value for F-statistic: < 2.2e-16\), który wskazuje, że model jest istotny statystycznie, że nie musimy wykluczyć zależności liniowej.
avPlots(model1)
Założenie 2. Zerowa średnia:
summary_model1_residuals <- skim(model1$residuals) %>% mutate(complete_rate= NULL) %>% mutate(n_missing= NULL) %>% mutate(numeric.hist = NULL) #niepotrzebne statystyki na ten moment
summary_model1_residuals
| Name | model1$residuals |
| Number of rows | 908 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | mean | sd | p0 | p25 | p50 | p75 | p100 |
|---|---|---|---|---|---|---|---|
| data | 0 | 9.02 | -38.05 | -5.07 | 0.83 | 5.67 | 26.24 |
Można powiedzieć, że średnia jest równa 0. Sprawdźmy to również za pomocą testu t-Studenta:
t.test(model1$residuals)
##
## One Sample t-test
##
## data: model1$residuals
## t = -1.3262e-15, df = 907, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.5874179 0.5874179
## sample estimates:
## mean of x
## -3.969378e-16
\(P-value=1\) oznacza, że nie ma wystarczających dowodów, aby odrzucić hipotezę zerową, co sugeruje brak istotnych statystycznie różnic. Zatem wektor losowy reszt w modelu ma średnią zero.
Założenie 3. Liniowa niezależność zmiennych objaśniających:
corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)
Możemy zobaczyć z macierzy korelacji, że np. zmienna \(Superplasticizer\) i \(Water\) są bardzo mocno ze sobą skorelowane, zatem są od siebie zależne. Założenie to odrzuca nasz model i musimy postąpić z korektą.
#Tworzenie drugiego modelu
Spróbujemy ulepszyć model1, aby całościowo nasz nowy model był lepszy oraz spełniał założenia. Skoro wystąpił problem, że dane zmienne objaśniające w modelu są od siebie zależne, należy im się przyjrzeć i poczynić odpowiednie kroki:
model1 <- lm(ConcreteStrength ~ ., data = concrete_data)
summary(model1)
##
## Call:
## lm(formula = ConcreteStrength ~ ., data = concrete_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.050 -5.071 0.829 5.674 26.235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.627526 26.526229 0.401 0.6888
## Cement 0.116777 0.008190 14.259 < 2e-16 ***
## BlastFurnaceSlag 0.097201 0.009827 9.892 < 2e-16 ***
## FlyAsh 0.075877 0.012028 6.309 4.41e-10 ***
## Water -0.212400 0.041257 -5.148 3.23e-07 ***
## Superplasticizer 0.244278 0.096365 2.535 0.0114 *
## CoarseAggregate 0.005622 0.009247 0.608 0.5434
## FineAggregate 0.007508 0.010535 0.713 0.4762
## Age 0.178619 0.006530 27.352 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.059 on 899 degrees of freedom
## Multiple R-squared: 0.6878, Adjusted R-squared: 0.685
## F-statistic: 247.6 on 8 and 899 DF, p-value: < 2.2e-16
Widzimy, że \(Pr(>|t|)\) dla zmiennych \(CoarseAggregate\) oraz \(FineAggregate\) \(>0.05\). Są potencjalnymi kandydatami do wyrzucenia z modelu.
corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)
Widzimy, że istnieje multikolinearności (dla zmiennej \(Water\) i \(Superplasticizer\) najbardziej). Zastosujemy użycie współczynnika VIF na naszym modelu. VIF jest współczynnikiem, który definiujemy dla poszczególnego współczynnika zmiennej objaśnianej. Powstaje jako iloraz wariancji współczynnika \(\beta_i\) w pełnym modelu do jego wariancji w modelu, który ma tylko \(x_i\) jako zmienną objaśniającą. Możemy go wyrazić wzorem \[ VIF(\beta_i) = \frac{1}{1 - R^2_{X_j|X_{-j}} }\] gdzie przez \(R^2_{X_j|X_{-j}}\) oznaczamy współczynnik \(R^2\) modelu regresji, gdzie zmienną objaśnianą jest \(X_j\), a objaśniającymi są wszystkie pozostałe zmienne niezależne z modelu pierwotnego. Można zauważyć z postaci tego współczynnika, że VIF jest zawsze większy od 1, a coraz większe wartości VIF wskazują na coraz większy stopień kolinearności. Praktyka pokazuje, że zaalarmowani powinniśmy być, gdy VIF danej zmiennej przekracza 5, a poważnie martwić, gdy przekracza 10. Zobaczmy jak to się prezentuje u nas:
car::vif(model1)
## Cement BlastFurnaceSlag FlyAsh Water
## 7.298968 7.647387 6.668511 7.325227
## Superplasticizer CoarseAggregate FineAggregate Age
## 3.106059 5.539066 6.719696 1.065311
Wnioski z wyniku:
Pomimo,że \(Cement\), \(BlastFurnaceSlag\) mają powyżej wartości 5, mają przeciętnie wysoką korelacje dodatnią ze zmienną \(ConcreteStrength\).
Pomimo,że \(Water\), \(Flyash\) mają powyżej wartości 5, mają przeciętnie wysoką korelacje ujemną ze zmienną \(ConcreteStrength\).
Dla zmiennych \(CoarseAggregate\) oraz \(FineAggregate\) również wyszły wartości powyżej 5. Z racji, że obie te zmienne najmniej pasowały do modelu, możemy spróbować stworzyć nowy model bez nich.
model2 <- lm(ConcreteStrength ~ . -CoarseAggregate-FineAggregate, data = concrete_data)
summary(model2)
##
## Call:
## lm(formula = ConcreteStrength ~ . - CoarseAggregate - FineAggregate,
## data = concrete_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.483 -5.032 0.858 5.614 26.049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.627889 4.136370 6.921 8.51e-12 ***
## Cement 0.111785 0.004109 27.202 < 2e-16 ***
## BlastFurnaceSlag 0.091174 0.004802 18.986 < 2e-16 ***
## FlyAsh 0.069215 0.007494 9.236 < 2e-16 ***
## Water -0.236652 0.020964 -11.288 < 2e-16 ***
## Superplasticizer 0.228643 0.087622 2.609 0.00922 **
## Age 0.178599 0.006525 27.372 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.052 on 901 degrees of freedom
## Multiple R-squared: 0.6876, Adjusted R-squared: 0.6855
## F-statistic: 330.6 on 6 and 901 DF, p-value: < 2.2e-16
Widzimy, że pomimo tego, że model się nie polepszył ale również się nie pogorszył znacząco, \(Adj. R^2\)=0.6876 vs. 0.6878.
Możemy zaobserwować, że wszystkie zmienne mają silne powiązanie z model z wyjątkiem \(Superplasticizer\). Wyeliminujemy ją oraz wtedy również multikorelacje z \(Water\).
Pozostałe wartości \(p-value\) są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Statystyka F również potwierdza istnienie związku pomiędzy zmienną \(ConcreteStrength\), a zmiennynmi objaśniającymi.
Stwórzmy kolejny model:
model3 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+Age, data = concrete_data)
summary(model3)
##
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh +
## Water + Age, data = concrete_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.524 -5.217 0.850 5.649 26.627
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.945116 3.611198 9.40 <2e-16 ***
## Cement 0.115895 0.003808 30.44 <2e-16 ***
## BlastFurnaceSlag 0.096952 0.004275 22.68 <2e-16 ***
## FlyAsh 0.080109 0.006243 12.83 <2e-16 ***
## Water -0.270327 0.016575 -16.31 <2e-16 ***
## Age 0.179726 0.006532 27.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.081 on 902 degrees of freedom
## Multiple R-squared: 0.6853, Adjusted R-squared: 0.6835
## F-statistic: 392.8 on 5 and 902 DF, p-value: < 2.2e-16
Również możemy zauważyć,że \(R^2\) nie zmienił się zbytnio w stosunku do poprzedniego modelu.
\(R^2\) na poziomie \(0.6853\), oznacza iż model wyjaśnia około 68,5% obserwacji wartości zmiennej \(ConcreteStrength\).
Wartości p-value są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Statystyka F również potwierdza istnienie związku pomiędzy zmienną \(ConcreteStrength\), a zmiennynmi objaśniającymi.
Zobaczmy współczynnik VIF:
car::vif(model3)
## Cement BlastFurnaceSlag FlyAsh Water
## 1.570321 1.440448 1.788143 1.176643
## Age
## 1.060627
Współczynnik ten jest poniżej 5 dla każdej zmiennej, zatem zmienne objaśniające nie są ze sobą skorelowane.
Zobaczmy, czy możemy ulepszyć powyższy model. Spójrzmy na histogramy naszych zmiennych:
par(mfrow = c(2,3))
names <- c("Cement","BlastFurnaceSlag","FlyAsh","Water","Age","ConcreteStrength")
for (i in names){
hist(concrete_data[[i]],
main = paste("Histogram - ", i),
xlab = "Value", ylab = "Frequecy")
}
Możemy potraktować zmienną \(Age\) logarytmem ze względu na jej wykładniczy charakter (zmienne \(FlyAsh\) oraz \(BlastFurnaceSlag\) mają wartości = 0).
model4 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = concrete_data)
summary(model4)
##
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh +
## Water + log(Age), data = concrete_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.6524 -4.1841 -0.1371 4.2187 21.4391
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.551137 2.443128 4.319 1.74e-05 ***
## Cement 0.119674 0.002559 46.772 < 2e-16 ***
## BlastFurnaceSlag 0.090130 0.002865 31.456 < 2e-16 ***
## FlyAsh 0.063087 0.004183 15.081 < 2e-16 ***
## Water -0.264682 0.011029 -23.999 < 2e-16 ***
## log(Age) 9.633018 0.182725 52.719 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.096 on 902 degrees of freedom
## Multiple R-squared: 0.8581, Adjusted R-squared: 0.8574
## F-statistic: 1091 on 5 and 902 DF, p-value: < 2.2e-16
\(R^2\) (oraz \(Adj. R^2\)) znacząco podwyższyły się w porównaniu do poprzedniego modelu (\(R^2\)=0,6853 vs 0.8581)
\(R^2\) na poziomie \(0.8581\), oznacza iż model wyjaśnia około 85,8% obserwacji wartości zmiennej \(ConcreteStrength\).
Wartości p-value są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Statystyka F również potwierdza istnienie związku pomiędzy zmienną \(ConcreteStrength\), a zmiennynmi objaśniającymi.
Przeprowadźmy jeszcze raz analizę założeń modelu:
Założenie 1. Liniowa zależność między zmienną objaśnianą, a objaśniającymi:
Możemy stwierdzić na podstawie \(Adj. R^2=0.8574\) oraz \(p-value for F-statistic: < 2.2e-16\), który wskazuje, że model jest istotny statystycznie, że nie musimy wykluczyć zależności liniowej.
avPlots(model4)
Założenie 2. Zerowa średnia:
summary_model4_residuals <- skim(model4$residuals) %>% mutate(complete_rate= NULL) %>% mutate(n_missing= NULL) %>% mutate(numeric.hist = NULL)
summary_model4_residuals
| Name | model4$residuals |
| Number of rows | 908 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | mean | sd | p0 | p25 | p50 | p75 | p100 |
|---|---|---|---|---|---|---|---|
| data | 0 | 6.08 | -16.65 | -4.18 | -0.14 | 4.22 | 21.44 |
Można powiedzieć, że średnia jest równa 0. Sprawdźmy to również za pomocą testu t-Studenta:
t.test(model4$residuals)
##
## One Sample t-test
##
## data: model4$residuals
## t = 2.9355e-15, df = 907, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3959577 0.3959577
## sample estimates:
## mean of x
## 5.922423e-16
\(P-value=1\) oznacza, że nie ma wystarczających dowodów, aby odrzucić hipotezę zerową, co sugeruje brak istotnych statystycznie różnic. Zatem wektor losowy reszt w modelu ma średnią zero.
Założenie 3. Liniowo niezależne zmienne objaśniające:
corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)
car::vif(model4)
## Cement BlastFurnaceSlag FlyAsh Water
## 1.573182 1.435734 1.781191 1.155967
## log(Age)
## 1.038097
Współczynnik ten jest poniżej 5 dla każdej zmiennej, zatem zmienne objaśniające nie są ze sobą skorelowane.
Założenie 4. Homoskedastyczność:
plot(model4)
bptest(model4)
##
## studentized Breusch-Pagan test
##
## data: model4
## BP = 10.529, df = 5, p-value = 0.06156
Otrzymaliśmy \(p-value\) znacznie większe niż 0.05, zatem nie mamy wystarczająco dowodów aby odrzucić hipotezę o homoskedastycznośći reszt.
Założenie 5. Normalność reszt:
Zbadamy to za pomocą histogramu reszt, wykresy kwantyl-kwantyl oraz testu Shapiro-Wilk’a:
ggplot(model4, aes(x=.resid)) + geom_histogram(bins=40) +labs(title = "Histogram reszt w modelu", x = "Wartość", y = "Częstotliwość")
ggplot(model4, aes(sample =.resid)) + geom_qq() + geom_qq_line() + labs(title = "Wykres kwantyl-kwantyl reszt", x = "Kwantyle teoretyczne", y = "Kwantyle próbkowe")
Histogram, oraz wykres kwantyl-kwantyl sugerują iż rozkład reszt może
być normalny. Zweryfikujemy to za pomocą testu Shapiro-Wilk’a:
shapiro.test(model4$residuals)
##
## Shapiro-Wilk normality test
##
## data: model4$residuals
## W = 0.99769, p-value = 0.237
Otrzymaliśmy \(p-value>0.05\), zatem nie mamy podstaw aby odrzucić hipotezę o normalności reszt.
Podsumowując, reszty w modelu spełniają klasyczne założenia modelu regresji.
Porównajmy nasz model ostateczny z poprzednimi:
st_err1 = summary(model1)$sigma
st_err2 = summary(model2)$sigma
st_err3 = summary(model3)$sigma
st_err4 = summary(model4)$sigma
paste("Residual Standard Error :", "Model 1:",round(st_err1,4),"|", "Model 2:",round(st_err2,4),"|", "Model 3:",round(st_err3,4), "Model 4:",round(st_err4,4) )
## [1] "Residual Standard Error : Model 1: 9.0591 | Model 2: 9.0516 | Model 3: 9.0807 Model 4: 6.0963"
Model o mniejszym błędzie standardowym jest lepszy od innych, czyli u nas model ostateczny - model4.
Dodatkowo jednym z czynników porównującym model może być również współczynnik determinacji \(Adj. R^2\), który tłumaczy jaki model lepiej wyjaśnia zmienność danych:
rsqr1 = summary(model1)$adj.r.squared
rsqr2 = summary(model2)$adj.r.squared
rsqr3 = summary(model3)$adj.r.squared
rsqr4 = summary(model4)$adj.r.squared
paste("Adjusted R^2 :", "Model 1:",round(rsqr1,4),"|", "Model 2:",round(rsqr2,4),"|", "Model 3:",round(rsqr3,4), "Model 4:",round(rsqr4,4) )
## [1] "Adjusted R^2 : Model 1: 0.685 | Model 2: 0.6855 | Model 3: 0.6835 Model 4: 0.8574"
Porównując współczynniki determinacji również możemy stwierdzić, że model4 najlepiej wyjaśnia zmienność danych.
Spróbujemy dodatko podzielić zbiór naszych danych na zbiór treningowy oraz zbiór testowy. Następnie porównamy nasz model z tym:
set.seed(1311)
split_data <- createDataPartition(concrete_data$ConcreteStrength, p = 0.75, list = FALSE)
train_data <- concrete_data[split_data, ]
test_data <- concrete_data[-split_data, ]
model4train <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = train_data)
summary(model4train)
##
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh +
## Water + log(Age), data = train_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.6395 -4.2858 -0.0804 4.3773 20.9613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.537399 2.790595 3.059 0.00231 **
## Cement 0.119369 0.002933 40.695 < 2e-16 ***
## BlastFurnaceSlag 0.089144 0.003391 26.285 < 2e-16 ***
## FlyAsh 0.063314 0.004836 13.093 < 2e-16 ***
## Water -0.251251 0.012638 -19.880 < 2e-16 ***
## log(Age) 9.519809 0.210946 45.129 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.138 on 677 degrees of freedom
## Multiple R-squared: 0.8553, Adjusted R-squared: 0.8543
## F-statistic: 800.5 on 5 and 677 DF, p-value: < 2.2e-16
Dodatkowym sposobem na poprawienie jakości modelu jest jeszcze spojrzeć na odległośći Cooka oraz dźwignie:
ggplot(model4train, aes(.hat, .stdresid))+ geom_point(aes(size=.cooksd))+ stat_smooth(method='loess', formula=y~x, se=FALSE)+ labs(title='Leverage vs Standardized Residuals',x='Leverage', y='Standardized Residuals', size='Cooks distance')
Postaramy się usunąć obserwacje o wysokiej dźwigni, oraz odległości
Cooka:
train_data2 <- train_data%>% mutate(cooks_D = cooks.distance(model4train))%>% filter(cooks_D <= 3*mean(cooks_D))
Sprawdzimy czy nasz model na zmienionych danych:
model4train2 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = train_data2)
summary(model4train2)
##
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh +
## Water + log(Age), data = train_data2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.8816 -3.8364 -0.1662 3.7222 15.4607
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.650197 2.580683 4.127 4.18e-05 ***
## Cement 0.118353 0.002617 45.231 < 2e-16 ***
## BlastFurnaceSlag 0.087462 0.003046 28.714 < 2e-16 ***
## FlyAsh 0.068838 0.004347 15.834 < 2e-16 ***
## Water -0.259741 0.011980 -21.680 < 2e-16 ***
## log(Age) 9.373795 0.187549 49.980 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.196 on 619 degrees of freedom
## Multiple R-squared: 0.8899, Adjusted R-squared: 0.889
## F-statistic: 1000 on 5 and 619 DF, p-value: < 2.2e-16
#Ostateczne porównanie modeli:
st_err4 = summary(model4)$sigma
st_err4train = summary(model4train2)$sigma
paste("Residual Standard Error :", "Model 4:",round(st_err4,4),"|", "Model 4 treningowy:",round(st_err4train,4))
## [1] "Residual Standard Error : Model 4: 6.0963 | Model 4 treningowy: 5.1959"
Model o mniejszym błędzie standardowym jest lepszy, czyli u nas model4train2.
rsqr4 = summary(model4)$adj.r.squared
rsqr4train = summary(model4train2)$adj.r.squared
paste("Adjusted R-squared :", "Model 4:",round(rsqr4,4),"|", "Model 4 treningowy:",round(rsqr4train,4))
## [1] "Adjusted R-squared : Model 4: 0.8574 | Model 4 treningowy: 0.889"
Również współczynnik determinancji jest wyższy w modelu model4train2, zatem uwzględnienie odległości Cooka polepszyło nasz model, oraz jego parametry na zbiorze testowym.
Najlepsze dopasowanie naszego modelu, przy jednoczesnym spełnieniu założeń regresji wielorakiej otrzymaliśmy w modelu ze zmiennymi objaśniającymi: \(Cement\), \(BlastFurnaceSlag\), \(FlyAsh\), \(Water\) oraz \(log(Age)\).
Dla tego samego modelu, poprzez nauczenie na zbiorze treningowym oraz uwzględnienie odległości Cooka otrzymaliśmy najlepsze dopasowanie, oraz najniższy bład standardowy na zbiorze testowym. Nasz model końcowy wyjaśnia około 88,9% zmiennej \(ConcreteStrength\) w zbiorze testowym, co możemy uznać za bardzo dobre dopasowanie.